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C/5 : ABSTRACT 
O 

The ideal gas equation of state with a constant adiabatic index, although 
[t] ! commonly used in relativistic hydrodynamics, is a poor approximation for most 

relativistic astrophysical flows. Here we propose a new general equation of state 
for a multi-component relativistic gas which is consistent with the Synge equa- 
tion of state for a relativistic perfect gas and is suitable for numerical (special) 
relativistic hydrodynamics. We also present a multidimensional relativistic hy- 
drodynamics code incorporating the proposed general equation of state, based 
on the HLL scheme, which does not make use of a full characteristic decomposi- 
tion of the relativistic hydrodynamic equations. The accuracy and robustness of 
this code is demonstrated in multidimensional calculations through several highly 
\ relativistic test problems taking into account nonvanishing tangential velocities. 

Results from three-dimensional simulations of relativistic jets show that the mor- 
[ phology and dynamics of the relativistic jets are significantly influenced by the 

different equation of state and by different compositions of relativistic perfect 
gases. Our new numerical code, combined with our proposed equation of state is 
very efficient and robust, and unlike previous codes, it gives very accurate results 
for thermodynamic variables in relativistic astrophysical flows. 
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Introduction 



Many high energy astrophysical phenomena, including accretion flows, jet flows, gamma- 
ray bursts, and pulsar winds involve relativistic flows. In powerful extragalactic radio sources, 
for example, ejections from galactic nuclei produce intrinsic beam Lorentz factors of usually 
more than 5 and apparently up to ~ 50, which are required to explain the apparent su- 
perluminal motions observed in extragalactic radio sources associated with active galactic 
nuclei (e.g.. lLister et al.ll2009l ). In the expansion of many relativistic jets the internal thermal 
energy of a gas is converted into bulk kinetic energy so as to reach a high Lorentz factor 
in a short distance. Then this kinetic energy is dissipated by shock interactions, mostly by 
terminal shock complexes, and partially by int ernal shocks within the jets as they propagate 



over long distances (e.g., iNorman et al.lll982l ). Since relativistic flows are inherently non- 



linear and complex, in addition to possessing large Lorentz factors, numerical simulations 
have been performed to investigate such relativistic flows, for example, in the propagation 
of relativistic extragalact i c jets (e.g., Duncan Hugheslll994l ; iMartf et al.lll997t iRosen et al. 



19991 : iHughes et al. 



2002: Choi et al. 



20071 ) 



Many explicit finite difference schemes originally applied to classical hydrodynamics 
have been employed to treat special relativistic hydrodynamics numerically. These schemes 



solutions to the local Riemann problem ( 


Schneider et al. 


19931: 


Falle & Komissarov 


1996|; 


Donat et al. 


1998 


Alov et al. 


1999; 


Del 


Zanna & Bucciantini 


2002 


Anninos & Fragile 


2003; 


Mienone & Bodo 


2005 




Mienone et al. 


2005 


). Some of these schemes adopt local charac- 



teristic decomposition of the Jacobian matrix to build numerical fluxes. It is often dif- 
ficult to build characteristic decomposition in some regimes, especially in ultrarelativistic 
limits, due to the degeneracy of the characteristic information. Thus the use of an alterna- 
tive scheme becomes sensible when the characteristic decomposition is unknown. The HLL 
scheme proposed by lHarten et al.l (119831 ) for classical hydrodynamics is based on an approx- 
imate Riemann solver that does not require full, and numerically expensive, characteristic 
decomposition. This feature of the HLL scheme makes its use very attractive, especially 
in multidimensions, where computational efficiency and robus tness is extreme l y imp ortant. 
This scheme was applied first to relat i vistic hydrodynamics by lSchneider et al.l (119931 ) in one 
dimension and by iDuncan fc Hughes! (119941 ) in multidimensions. 



It is worth stressing that many treatments of relativistic astrophysical problems have 
assumed a ideal gas equation of state with a constant polytropic index, but this is a reasonable 
approximation only if the gas is either strictly subrelativistic or ultrarelativistic. However, 
when the gas is semirelativistic or when the gas has two components, e.g., nonrelativistic 
protons and relativistic electrons, this assumption is no longer correct. This was shown for 



-3 - 



the relativistic perfect gas law by ISyngd ( 119571 ). where the exact form of an equation of 
state relating thermodynamic quantities of specific enthalpy and temperature is completely 
described in terms of modified Bessel functions. 

Since the correct equation of state for the relativistic perfect gas has been recognized 
as being important, several investigations with a more general equatio n of st ate have been 
reported in numerical relativistic hydrodynamics. iFalle fc Komissarovl (119961 ) described an 
upwind numerical code for special relativistic hydrodynam ics with the Synge eq uation of 
state for multi-component relativistic gas. More recently, iMignone et al.l (120051 ) used, in 
their upwind relativistic hydrodynamic code, a simple equation of state that closely ap- 
proximates the Synge equation of state for a single-component relativistic gas. Several 
numerical simulations in the context of relativistic extragalactic jets make use of the gen- 
eral eq uation of state to account for transitions from nonrelativistic to relat i vistic temper- 



ature (Komissarov fc Falle 



1998 



Scheck et al. 2002: Perucho fc Marti 2007: Meliani et al. 



20081 : iRossi et al.l 120081 ) . IScheck et al.l (120021 ) used the Synge equation of state for differ- 
ent compositions, including pure leptonic and baryonic plasmas, to investigate the influence 
of the compos i tion o f relativistic extragalactic jets on their long-term evolution. Similarly, 
Meliani et al.l (120081 ) studied the relativistic extragalactic jet deceleration through density 
discontinuities by using the Synge-like equation of state with a variable polytropic index. 

In this work we propose an analytical form of equation of state for the multi-component 
relativistic perfect gas that is consistent with the Synge equation of state in the relativistic 
regime. This proposed equation of state is suitable for a numerical code from the computa- 
tional point of view, unlike the Synge equation of state, which involves the computation of 
Bessel functions. We then build a multidimensional relativistic hydrodynamics code based 
on the HLL scheme using the proposed equation of state, and demonstrate the accuracy and 
robustness of this code by presenting several test problems and numerical simulations. In 
particular, we plan to use this code to simulate relativistic extragalactic jets that are proba- 
bly composed of a mixture of relativistic particles of different masses. Numerical simulations 
of relativistic jets of different compositions are challenging, but are made tractable by using 
the proposed equation of state that accounts for different compositions of relativistic gas. 

This paper is organized as follows. In §2 we present the equations of motion and the 
Synge equation of state, and we describe the proposed general equation of state for the 
relativistic gas. In §3 we describe a relativistic hydrodynamics code based on the HLL 
scheme, incorporating this proposed equation of state. In §§4 and 5 we present numerical 
tests and simulations with the code to demonstrate the performance of the code. A conclusion 
is given in §6. 
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2. Relativistic Hydrodynamics 

2.1. Equations of Motion 

The motion of relativistic gas is described by a system of conservation equations. The 
equations in special relativistic h ydrodynamics are written in a covariant form (ILandau &: Lifshitz 
19591 : Iwilson fc Mathewsl 12003 1) as 

d a (pU a ) = 0, (1) 

d a (phU a U p + pg aP ) = 0. (2) 

Here, d a = d/dx a is the covariant derivative with respect to spacetime coordinates x a = 
[t,Xj], U a = [T,Tvj] is the normalized {U a U a = — 1) four- velocity vector, where T is the 
Lorentz factor, and g af5 = diag{ — 1, 1, 1, 1} is the metric tensor in Minkowski space. The 
speed of light is set to unity (c = 1) in this work. Greek indices (e.g., a, (3) denote the 
spacetime components while Latin indices (e.g., i, j) indicate the spatial components. The 
rest mass density, specific enthalpy, and pressure in the local rest frame are denoted by p, h, 
and p, respectively. 

In Cartesian coordinates these relativistic hydrodynamic equations can be written in 
conservative form as 

8q . dF, dF„ dF; 



+ 



+ 



0. 



dx dy dz 
where q is the state vector of conservative variables and F 
the flux vectors in the x, y, and ^-directions, defined by 



F 

xi 1 yi 



(3) 

and F z are respectively 



" D - 




Dv x 


M x 




M x v x + p 


My 


F = 


M y v x 


M z 




M z v x 


. E . 




_(E + p) v x _ 



(4) 



The flux vectors F y and F z are given by properly permuting indices. The conservative 
variables D, M x , M y , M z , and E represent respectively the mass density, three components 
of momentum density, and energy density in the reference frame. 

The variables in the reference frame are nonlinearly coupled to those in the local rest 
frame via the transformations 

D = Tp, (5) 



M x = T phv x , M y = T phvy, 
E = T 2 ph - p, 



M z = T 2 phv z 



(6) 
(7) 



where the Lorentz factor is given by T = 1/Vl - v 2 with 



v 2 + v 2 + v 2 . 

V x I Uy I U z . 
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2.2. Equation of State 

The system of conservation equations describing the motion of relativistic gas is com- 
pleted with an equation of state that relates the thermodynamic quantities of specific en- 
thalpy, rest mass density, and pressure. In general the equation of state can be expressed 
with the specific enthalpy expressed as a function of the rest mass density and the pressure 

h = h(p,p). (8) 

The sound speed c s is then defined as 

2 pdh ( dh\ _1 



The explicit form of the sound speed depends on the particular choice of the equation of 
state. 

The exact form of equatio n of state for a relativistic perfect gas composed of multi- 



ple components was derived by ISyngd (119571 ). For a single-component relativistic gas the 
equation of state is written as 

h = iwr (10) 

where K 2 and K 3 are respectively the modified Bessel function of the orders two and three 
and £ = p/p is a measure of inverse temperature. Under the Synge equation of state a 
relativistic perfect gas is entirely described in terms of the modified Bessel functions. The 
Synge equation of state for a relativistic gas can be written in the form 

h = l + (11) 

X - 1 p 

where the quantity 7* is defined by 

7* = k ~ l • (12) 

Then the sound speed is written as 

2 IrP M „x 

c * = w (13) 

and the relativistic adiabatic index 7 r is given by 

h'£ 2 

i' = WTi> (14) 

where h' = dh/d^. The quantities 7* and 7 r are constant and equal if the gas remains 
ultrarelativistic or subrelativistic (i.e., 7* = 7 r = 4/3 for £ 1 or 5/3 for £ ^> 1). For 
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the intermediate regime 7 * and 7 r vary slightly differently between the two limiting cases 
( IFalle fc Komissarovlll996l ). 



For a multi-component relativistic gas, the direct use of the Synge equation of state 
involves the computation of Bessel functions, thus requires significant computation cost and 
results in computational inefficiency Here we propose a new general equation of state for 
mult i- component relativistic gas that uses analytical expression and is more efficient and 
suitable for numerical computations. We suppose that the relativistic gas is composed of 
electrons, positrons, and protons although more components of relativistic gas easily can be 
considered. The total number density n is then given by 



n 



71: 



+ n e + + n. 



(15) 



where n e -, n e +, and n p + are the electron, positron, and proton number densities, respec- 
tively. We ignore the production or annihilation of electron-positron pairs and assume the 
composition of electrons, positrons, and protons is maintained. The assumption of charge 
neutrality gives us the relations n e - = n e + + n p + and n = 2n e ~ . The total rest mass density 
and pressure are respectively given by p = ^2 Pi = Y n irrii and p = YVi = Y nikT, where 
k is the Boltzmann constant and T is the temperature. 

For our equation of state for multi-componen t relativistic gas w e adopt the equatio n 
of state, previously introduced by iMathewsl (119711 ) and later used by iMeliani et al.l ( 12004 ). 
which closely reprodu ces the Synge equation of state for a single-component relativistic gas 
(IMignone et al.l 120051 ) . The equation of state takes the form 



V 
P 



(16) 




where e is the energy density in the local rest frame. It can be solved for the specific enthalpy 
using ph = e + p as 

c „ In /^\ 2 

(17) 

For a multi-component relativistic gas, the total enthalpy is then given by ph = ^2 pihi = 
^2[(5/2)riikT + a/ (9/A)(riikT) 2 + (r^m,) 2 ]. After we express the total enthalpy in terms of 
each component, we can eliminate n e + using the above relations given in the charge neutrality 
assumption. Our proposed equation of state is obtained by simplifying the expression of the 
total enthalpy. We find the resulting equation of state for a mult i- component relativistic gas 
to be 



h 



51 
2^ 



X) 



9 1 



+ 



1 



nl/2 



16 e (2-x + XPY 



+ x 



9 1 



+ 



16 e (2-x + XPY 



1/2 



-7- 



where \ = n p + / n e - is the relative fraction of proton and electron number densities and 
H = m p /m e is the mass ratio of proton to electron. We note that \ = represents an 
electron-positron gas while % — 1 indicates an electron-proton gas. The proposed equation of 
state reduces to the equation of state in equation (17) when the electron-positron gas (x = 0) 
is considered. This proposed equation of state holds only in the limit where the composition 
of the plasma is fixed in space and time. In reality, however, the plasma composition may 
change through fluid mixing or the creation or annihilation of electron-positron pairs, either 
of which can alter the equation of state by effectively changing \ in different regions of the 
fluid. Nonetheless, this assumption of constant \ is a useful first order approximation in 
many situations. 

Using the proposed equation of state we show in Figure [TJ the relativistic adiabatic 
indices, 7 r , as functions of inverse temperature £, for several different compositions of the 
relativistic gas. Compositions of x — 0, 0.1, 0.3, 0.6, and 1 are shown. The value of 7 r 
critically depends on the composition of the relativistic gas as protons become relativistic at 
a much higher temperature than do electrons. Different proton fractions cause the adiabatic 
index to vary substantially at intermediate temperatures. For electron-positron and electron- 
proton gases, we compare 7 r values from our equation of state with 7 values from the Synge 



equation of state as given in Figure 1 of iFalle fc Komissarovl ( 119961 ). Direct measurements 



of relative errors give \j r — 7I/7 < 0.4% for an electron-positron gas and < 0.3% for an 
electron-proton gas. 

Figure |5] shows the quantities 7*, specific enthalpy h, and sound speed c s , as functions 
of inverse temperature £, for different equations of state. Results are shown for the ideal gas 
equation of state with fixed adiabatic indices of 7 = 5/3 and 4/3 as well as for the proposed 
general equation of state for pure electron-positron and pure electron-proton gases. The 
thermodynamic quantities computed using the proposed general equation of state asymp- 
totically approach those computed using the ideal gas equation of state with 7 = 4/3 in 
the hot gas limit < 1) and 7 = 5/3 in the cold gas limit (£ ^> 1). For intermediate 
regimes the thermodynamic quantities vary between those two limiting cases, depending on 
the composition of the relativistic gas. 



Numerical Scheme 



We now describe a relativistic hydrodynamics code incorporating the proposed equa- 
tion of state f or a r elativistic perfect gas, based on the HLL scheme originally proposed by 
Harten et al.l (119831 ) for classical hydrodynamics. The HLL scheme avoids a full characteristic 
decomposition of the relativistic hydrodynamic equations and uses an approximate solution 
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to the Riemann problem where the two constant states separated by the middle contact wave 
are averaged into a single intermediate state. The HLL scheme requires accurate estimates 
of the maximum and mi nimum wave sp eeds for the solution of the Riemann problem. In 



classical hydrodynamics lEinfeldtl ( 119881 ) proposed good ways to estimate the wave speeds 



based on the maximum and minimum eigenvalues of the Jacobian matrix. 

The wave speeds needed in our formulation for relativistic flows can be similarly es- 
timated from the maximum and minimum eigenvalues of the Jacobian matrix, A x = 
dF x (u)/dq(u), where u = [p, v x , v y , v z ,p] T is the state vector of primitive variables, and 



;i - eg) v x ± V(i-« 2 )c* [i-v*c* -(i-<* 



21 



1 — v z c z s 

The maximum and minimum eigenvalues are based on a simple application of the relativistic 
addition of velocity components decomposed into coordinates directions and simply reduce 
to a x = (v x ± c s )/(l ± v x c s ) in the one- dimensional case. 

Numerical integration of relativistic hydrodynamic equations advances by evolving the 
state vector of conservative variables in time. However, in order to compute the flux vectors 
for the evolution, the primitive variables u involved in the flux vectors should be recovered 
from the conservative variables q at each time step by the inverse transformation 

D , s 

P = p, (20) 

M x M y M z 

v x = t; 1 v v = — — — , v z = — , (21) 

E + p' y E + p' E + p 1 v ; 

p = YDh - E. (22) 

The inverse transformation is nonlinearly coupled and reduces to a single equation for the 
pressure 

f(p) = T(p)Dh(ap))-E-p = 0. (23) 

This nonlinear equation can be solved numerically using a Newton- Raphson iterative method 
in which the derivative of the equation with respect to pressure is given by 

df dT ^jjdhd£ ^ 
dp dp dt; dp 

Once the pressure is found numerically, the rest mass density and velocity are recovered 
by the inverse transformation. This procedure of inversion from conservative to primitive 
variables is valid for a general equation of state by specifying the expression of specific 
enthalpy. 
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The numerical integration of relativistic hydrodynamic equations proceeds on spatially 
discrete numerical cells in time, based on the finite difference method. In our implementation 
the state vector q n at the cell center i at the time step n is updated by calculating the 
numerical flux vector f™*+y 2 along the re-direction at the cell interface i + 1/2 at the half 
time step n + 1/2 as follows 



At r 

„n+l _ „n ^ L 



Ax 



( -n+1/2 _ -n+1/2 \ (C) ^ 

yj x,i+i/2 J x,i-i/2j ■ y* ) 



The numerical flux vector is calculated from the approximate Riemann solution and is given 
in the form 

a + F (u n+1/2 )-a~ F (u n+1/2 ) 

-n+1/2 a x,i+l/2 r x\ u L,i+l/2) a x,i+l/2 r *\ u R,i+l/2) 

J .r.i ■ I 2 ~ 



a x,i+l/2 a x,i+l/2 



a x,i+l/2 a x,i+l/2 



, n+1/2 x / n+1/2 \ 

q( u L,i+i/2) - 9(^+1/2) 



a x,i+l/2 ~ a x,i+l/2 

where the maximum and minimum wave speeds are defined by 

,+ _ ^ r n + /„ n+1/2 \ +/ n+1/2 



(26) 



"x.i+l/2 = maX { ' a t( U L,i+l/2)> a x( U R,i+l/2)}> ( 27 ) 

+1/2 = min{0, a-(u n L ^ /2 ), a~ (u£^ /2 )} . (28) 



a 

X,l 



Here w2t+i/2 an< ^ u r + h-i/2 are ^ ne an< ^ right state vectors of the primitive variables, 
which are defined at the left and right edges of the cell interface, respectively. In the first 
order of spatial accuracy the left and right state vectors reduce to w^+1/2 = u i +1 ^ 2 an d 
u R + i+i/2 = U i+i^- For second-order accuracy in space the left and right state vectors are 
interpolated as 

_ n+1/2 1 A n n+1/2 _ n+1/2 _ 1 « „ f2q x 

with the minmod limiter 

A< = 1 [sign(A + <) + sign(A-<)] min{|A+<| , |A"<|}, (30) 

or the monotonized central limiter 

A< = l - [sign(A+<) + sign(A-<)] min{2 | A+<| , 2 |A~<| , l - |A+< + A~<|}, (31) 

where A + w™ = u™ +1 — and A~u n = uf — uf^. In our formulation, the state vector of the 
primitive variables defined at the half time step, ■u™ +1//2 , is computed from a predictor step, 
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q n+l/2 = € _ (Ar/2Ax)(/ n . +i/2 _ / » i _ i/2 ) > with the flux vector f n, +i/2 calculated by 

replacing the time step n + 1/2 by the time step n in equations (26) to (29). This approach 
makes our code second order in time as well. 

The time step is restricted by the Courant condition At n = C cour Ax/max{|a^ +1 , 2 |} 
with C cour < 1. For multidimensional extensions, the numerical integration along the re- 
direction is applied separately t o the y- (and z-) directions through the Strang-type di- 
mensional splitting ( IStrana Il968l ) . In order to maintain second-order accuracy the order of 
dimensional splitting is completely permuted in each successive sequence. 



4. Test Problems 

We have applied the numerical scheme with the general equation of state, described in 
previous sections, to several test problems. In all the test problems the minmod limiter and 
a Courant constant, C cour = 0.4, are used. 



4.1. Relativistic Shock Tube 



The relativistic shock tube test is characterized initially by two different states separated 
by a discontinuity. As the initial discontinuity decays, distinct wave patterns consisting of 
shock waves, contact discontinuities, and rarefraction waves appear in the subsequent flow 
evolution. In the relativistic shock tube problem the decay of the initial discontinuity signifi- 
cantly depends on the tangential velocity since the velocity components are coupled through 
the Lore ntz factor in the e quations and the specific enthalpy also couples with the tangential 
velocity (IPons et al.ll2000l ). As a result, the relativistic shock tube problem becomes more 
challenging in the presence of tangential velocitie s. This relativist i c shock tube problem is 
very good test since it has an analytic solution (IPons et al.l |2000| ; iGiacomazzo fc Rezzolla 
20061 ) where the numerical solution can be compared. 



We have performed two sets of the two-dimensional relativistic shock tube tests using the 
general equation of state for an electron-positron gas. The first test set is less relativistic, 
with h ~ 5, while the second test set is more severe, with a large initial internal energy 
(h ^> 1). For each of the test sets, two cases are presented, one having only parallel velocity 
components while the other has, in addition, tangential velocities. In the first set, the initial 
left and right states for the case of only parallel velocities are given by pi = 10, pn = 1, 
v p ,l = 0, v Pj r = 0, pi = 13.3, and pr = 10~ 6 ; for the case where tangential velocities are 
included these additional initial values are v t l — 0.9 and v t R = 0.9. For the second set 
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the initial left and right states for the case of only parallel velocity are pi — 1, pr — 1, 
v p,l = 0, v p> r = 0, pl = 10 3 , and pr = 10~ 2 ; now for the case when tangential velocities are 
considered as well, v^l = 0.99 and v t) R = 0.99 are taken. Here the subscripts L and R denote 
the left and right states separated by an initial discontinuity placed along the main diagonal 
in the two-dimensional computational plane. Structures such as waves and discontinuities 
propagate along the diagonal normal to the initial discontinuity. Here v p is velocity parallel 
to the wave normal in the plane, given by v p = ^Jv 2 + v 2 y and v t is velocity tangential to the 
wave normal in the direction out of the plane, given by Vt = v z . The numerical computations 
are performed in a two-dimensional box with x = [0, 1] and y = [0, 1] using 512 x 512 cells 
for the first set and 2048 x 2048 cells for the second set. 

The results from the numerical computations carried with the general equation of state 
for an electron-positron gas for the first test set are shown in Figure [21 Wave structures 
are measured along the main diagonal line x = y = to 1 at times t = 0.4a/2 and 0.8a/2 
for the cases of the only parallel and also tangential velocities, respectively. The numerical 
solutions are marked with open circles and the ana l ytical solutions obtained using the nu- 



merical code available from iGiacomazzo fc Rezzollal (120061 ) are plotted with solid lines. Our 



numerical scheme with the proposed general equation of state is able to reproduce all the 
wave structures with very good accuracy and stability, as shown in Figure El The shock 
waves and rarefraction waves are captured correctly, while the contact discontinuities are 
relatively more smeared due to the use of the HLL scheme and the minmod limiter. The 
inclusion of the tangential velocity leads to the same basic wave pattern as in the absence of 
tangential velocity, but the numerical solutions are significantly modified. 

Figure H] shows the results from the numerical computations with the general equation of 
state for an electron-positron gas for the second test set. Structures are measured along the 
main diagonal line x — y — to 1 at times t = 0.4a/2 and 1.8v^2 for the cases of the parallel 
only and included tangential velocities, respectively. The numerical solutions compare well 
with the analytical solutions. Shock waves and contact discontinuities propagate to the right, 
while rarefraction waves move to the left. As shown in Figure HI all the wave structures 
are accurately reproduced and their stability is good; however, the contact discontinuities 
are somewhat smeared. Again, the inclusion of the tangential velocity has a considerable 
influence on the numerical solutions. 

For a quantitative comparison with analytical solutions we have calculated the norm 
errors of the rest mass density, parallel and tangential velocities, and pressure defined by, 
e.g., for density, = ^ . — p^\AxiAyj, where the superscripts N and A represent 

numerical and analytical solution, respectively. The norm errors are given in Table [H and 
demonstrate a very good agreement between the numerical and analytical solutions in all 
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the primitive variables. We also carried out the two-dimensional relativistic shock tube tests 
considered in Figures [3] and H]for several other possible combinations of the tangential velocity 
pairs v tj L = 0, 0.9, 0.99 and v tj R = 0, 0.9, 0.99. In general, the stable wave structures are 
reproduced and the numerical solutions are reasonably comparable to the analytical solutions 
for all the different combinations of tangential velocity pairs. As a consequence, the results 
from the two-dimensional relativistic shock tube tests show that our code is able to robustly 
and accurately capture discontinuities and waves. 



4.2. Relativistic Shock Reflection 

The relativistic shock reflection problem involves a collision between two equal gas flows 
moving at relativistic speeds in opposite directions. The collision of the two cold gases causes 
compression and heating of the gases as kinetic energy is converted into internal energy. This 
generates the two strong shock waves to propagate in opposite directions, leaving the gas 
behind the shocks stat ionary. The analytical sol ution of this relativistic shoc k refle ction 



problem is obtained by iBlandford fc McKed (119761 ) and iGiacomazzo fc Rezzollal (120061 ) 



We have tested the two-dimensional relativistic shock reflection problem using the gen- 
eral equation of state for an electron-positron gas. Two cases are presented here. As for the 
shock tube tests, one includes only the parallel velocity and the other includes in addition 
the tangential velocity. The left and right states for the case of only parallel velocity are 
initially given by pl — 1, Pr — 1, v Pt L = 0.99, v Pt R = —0.99, pl = 10~ 6 , and pr = 10~ 6 , and 
for the case when the tangential velocity is included, oppositely directed tangential velocities 
v t ,L = 0.1 and v t> R = —0.1 are assumed to be present on either side of the plane. Here the 
subscripts L and R stand for the left and right states separated by the initial collision points 
located along the main diagonal in the two-dimensional computational plane. The two shock 
waves propagate diagonally in opposite directions, and should keep symmetric with respect 
to the initial collision points. As in the relativistic shock tube tests, v p is velocity parallel to 
the wave normal in the plane and v t is velocity tangential to the wave normal in the direction 
out of the plane. The numerical computations are carried out in a two-dimensional box with 
x = [0, 1] and y = [0, 1] using 512 x 512 cells. 

Figure [5] shows the results from our relativistic shock reflection tests with the general 
equation of state for an electron-positron gas. Structures are measured along the main 
diagonal line x — y — to 1 at time t = 0.8a/2. The numerical solutions are in very 
good agreement with the analytical solutions. In both cases, the shock wave is resolved by 
two numerical cells, and there are no numerical oscillations behind the shocks. As shown 
in Figure the compression ratio between shocked and unshocked gases is about 30 for 
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the rest mass density and about 70 for the pressure in the case of only parallel velocities; 
the inclusion of the tangential velocities increases the compression ratio to about 40 for 
the rest mass density and to about 140 for the pressure. Near x = y = 0.5, the density 
distribution slightly underestimates the analytical solution and a stationary discontinuity in 
the tangential velocity is somewhat diffused due to the numerical effect of reflection heating 
phenomena. The norm errors of the rest mass density, parallel and tangential velocities, and 
pressure are also given in Table [H A direct comparison with the analytical solutions shows 
that the measured errors are very small for all the primitive variables. 

The accuracy of numerical solutions depends on the number of cells spanned by the 
computational box. We have run the relativistic shock reflection test in Figure [5]^b) with 
different numerical resolutions to check the convergence rate. Except for the numerical 
resolutions the initial conditions are identical to those used in the test in Figure E(b). We 
have computed the norm errors for rest mass density, velocities, and pressure with different 
resolutions. Numerical resolutions of 16 2 , 32 2 , 64 2 , 128 2 , 256 2 , and 512 2 cells give norm errors 
for the rest mass density of 7.56, 4.74, 2.38, 1.17, 0.48, and 0.25, respectively. As expected 
for discontinuous problems, first-order convergence in the norm errors for rest mass density is 
obtained with increasing the numerical resolution. Similar clear trends toward convergence 
are seen in the norm errors for velocities and pressure. 



4.3. Relativistic Blast Wave 

In a relativistic blast wave test a large amount of energy is initially deposited in a 
small finite spherical volume and the subsequent expansion of that overpressured region is 
evolved forward in time. This produces a spherical shock propagating outward from an initial 
discontinuity at an arbitrary radius. This radial blast wave explosion provides a useful test 
problem to explore the spherically symmetric properties in highly relativistic flow speeds. 

We have performed the three-dimensional relativistic blast wave test with the general 
equation of state for an electron-proton gas. The initial condition for the relativistic blast 
wave problem consists of two constant states given by pi — 1, po = 1, vi = 0, vo = 0, 
pi = 10 3 , and po = 1, where subscripts / and O represent the inner and outer states separated 
by an initial discontinuity at the radius r = 0.5 in the three-dimensional computational box. 
The numerical computations are performed in a three-dimensional box with x = [0,1], 
y = [0,1], and z = [0,1] using 256 x 256 x 256 cells. Outflow boundary conditions are 
used at all boundaries except at the symmetry axis where reflecting boundary conditions are 
imposed. 
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Figure [6] shows the profiles of rest mass density, radial velocity, and pressure along the 
radial distance in the relativistic blast wave tests. Radial structures are measured along the 
main diagonal line x = y = z = 0tola.t time t = 0.4. Numerical computations carried 
out with our general equation of state for an electron-proton gas, as well as with the ideal 
gas equation of state with constant adiabatic index 7 = 5/3 are shown. In both cases, 
high Lorentz factors (up to ~ 15) are generated and the radial symmetry is well preserved. 
However, compared to the 7 = 5/3 case, the results obtained with the electron- proton case 
show significant differences. In fact, the spherical shock wave has a higher density peak and 
propagates at a slower radial velocity in the electron-proton case. 

In Figured we present the images of rest mass density, pressure, and Lorentz factor in 
the relativistic blast wave test with the general equation of state for an electron-proton gas. 
The images are shown in logarithmic scales on the plane x = at time t = 0.4. As shown 
in the images, the spherical shock wave propagates to larger radius very well, preserving the 
initial spherical symmetry. 



5. Relativistic Axisymmetric Jets 



As a practical astrophysical problem using this code, we consider the propagation 
of relativistic axisymmetric jets in three dimensions. The simulation of a relativistic jet 
has been presented as a test simulation in almost all relativistic hydro dy namic codes, us 



ing the ideal gas equation of state with a constant ad iabatic index (e.g., lAloy et al.lll999 



Del Zanna &: Bucciantinill2002l ; iMignone &: Bodoll2005l ). We have performed this relativistic 
jet simulation to confirm the accuracy and robustness of our numerical scheme incorporating 
the proposed general equation of state as well as to make a preliminary investigation of the 
influence of the equation of state on the propagation of relativistic jets. 

Some previous investigations of relativistic jet propagation ha ye concentrated on the 
i mpor tance of a realistic equation of state in relativistic jet flows. IMignone &: McKinney 



(120071 ) addressed the effect of a realistic equation of state on relativistic (magnetized) flows 
including relativistic jets. They showed that the choice of a realistic equation of state can 
significantly alter the solut ion when large temperat ure gradients are present. To study the 
evolution of low-power jets iPerucho fc MartJ (120071 ) performed numerical simulations using 
an equation of st ate for a two - compo nent relativistic gas that separately treats leptonic and 
baryonic matter. iRossi et all (120081 ) examined the dynamical evolution of relativistic light 
jets in the presence of an induced perturbation using a Synge-like equation of state for a 
single-component relativistic gas. 
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The approximate propagation velocity of the jet through the homogeneous ambient 
medium can be derived from momentum flux balance between the beam and ambient medium 
in th e reference frame of the working surface that separates beam and shocked ambient gas 
(e.g., iMartf et al.lll997l ). Assuming pressure equilibrium between the beam and the ambient 
medium, the one- dimensional jet advance velocity, estimated in the rest frame of the ambient 
medium, is then 

T b y/-qh h /h a 

Va = \ / I n Vb > ( 32 ) 

1 + LbVV h b/h a 



where rj and T b are given by rj = Pb/p a and T b = 1/ a/1 — v\. Here the subscripts b and a 
indicate the beam and the ambient medium, respectively. The morphology and dynamics of 
the relativistic jet propagating into the homogeneous medium is commonly specified by the 
beam to ambient medium density ratio rj, the beam Lorentz factor r&, and the beam Mach 
number M& = Vb/c St b- 

We have considered the three-dimensional simulations of relativistic jets using different 
equations of state. Table [2] lists the simulation parameters of the four different cases corre- 
sponding to the ideal gas equation of state with constant adiabatic index 7 = 5/3 and 4/3 
and the proposed general equation of state for electron-positron and electron-proton gases. 
In all these cases, we choose the beam density pb = 0.1, the ambient medium density p a = 10, 
the beam velocity Vb = 0.99 and we assume that the beam is in pressure equilibrium with 
the ambient medium, pb = p a . This gives the density contrast r] = 10 -2 , the beam Lorentz 
factor = 7.1, and the beam Mach number Mb = 2. The numerical simulations have 
been performed in the three-dimensional computational box with x = [0, 1], y = [0, 1], and 
z = [0,4] using a uniform numerical grid of 128 x 128 x 512 cells. The beam has an initial 
radius r b = 1/8 (corresponding to 16 cells), is launched from the origin, and propagates 
through a uniform static ambient medium along the positive z-direction. Outflow boundary 
conditions are set at all boundaries except along the symmetry axis where reflecting bound- 
ary conditions are used and in the injection region where an inflow boundary condition is 
imposed to keep the beam constantly fed. The monotonized central limiter and the Courant 
constant C COUT = 0.2 are used in all these jet propagation cases. 

The images in Figure [8] display the logarithms of the rest mass density on the plane x = 
at time t = 6 for the four different cases in the simulations of relativistic axisymmetric jets. 
The bow shock, the beam, and the cocoon surrounding the beam can be clearly identified in 
all four cases, confirming the ability of our code to follow complex relativistic flows. For each 
case, a bow shock that separates the shocked jet material from the shocked ambient medium 
is driven into the ambient medium. The beam itself is terminated by a Mach disk where 
much of the beam's kinetic energy is converted into internal energy. Shocked jet material 
flows backward along the working surface into a cocoon, resulting in the development and 
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mixture of turbulent vortices in the cocoon, and the interaction of these vortices with the 
beam forms oblique internal shocks within the beam close to the Mach disk, which causes 
the deceleration of the jet. The four different cases, however, show differences in specific 
morphological and dynamical properties. The cold jet (case A) propagates at the slowest 
velocity, and the jet produces a broad bow shock, a thick cocoon, and has the Mach disk 
located quite far behind the bow shock. On the contrary, the hot jet (case B) is dominated 
by a narrow bow shock, has a thin cocoon, and its Mach disk lies very close to the bow 
shock, thanks to its having the fastest advance velocity. The electron-positron and the 
electron-proton jets (cases C and D) propagate faster than the cold jet, but slower than the 
hot jet. Therefore, the electron-positron and the electron-proton jets possess morphological 
and dynamical properties intermediate between the cold and the hot jets. In terms of our 
simulation parameters, the electron-positron jet tends to be more similar to the cold jet, 
while the electron-proton jet seems to share more features with the hot jet. 

In Figure [9j the position of the bow shock is plotted as a function of time for the four 
different gases. The symbols mark the numerical estimate for a selected time interval for each 
case, and the lines represent the one- dimensional theoretical estimate in equation (32). The 
numerical simulations are in good agreement with the one- dimensional theoretical estimates 
for the bow shock location for all cases. 



6. Conclusion 

Most numerical codes for special relativistic hydrodynamics have used the ideal gas 
equation of state with a constant adiabatic index, but this is a poor approximation for most 
relativistic astrophysical flows. We proposed a new general equation of state for a multi- 
component r elativistic gas, based on the Synge equation of state for a relativistic perfect gas 



( ISyngd 119571 ). Our proposed general equation of state is very efficient and suitable for nu- 
merical relativistic hydrodynamics since it has an analytic expression. The thermodynamic 
quantities computed using the proposed general equation of state behave correctly asymp- 
totically in the limits of hot and cold gases. For intermediate regimes the thermodynamic 
quantities vary between those for the two limiting cases, depending on the composition of 
the relativistic gas. 

We also presented a multidimensional relativistic hydrodynamics code incorporating 
the proposed general equation of sta te for a multi-compo nent relativistic gas. Our numerical 



code is based on the HLL scheme (IHarten et al.lll983l ). which avoids a full characteristic 



decomposition of the relativistic hydrodynamic equations and uses an approximate solution 
to the Riemann problem for flux calculations. Since the numerical code is fully explicit, re- 
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taining a second- order accuracy in space and time, it is simple to extend the code to different 
geometries or to produce parallelized versions of this code. The analytical formulation of the 
proposed equation of state and the numerical scheme being free of complete characteristic 
decomposition make the code very efficient and robust in ultrarelativistic multidimensional 
problems. 

The accuracy and robustness of the code are demonstrated in two dimensions using the 
test problems of the relativistic shock tube and the relativistic shock reflection and in three 
dimensions using the test problem of the relativistic blast wave. The direct comparisons of 
numerical results with analytical solutions show that shocks and discontinuities are correctly 
resolved even in highly relativistic test problems with nonvanishing tangential velocities. 
Results from the three-dimensional simulations of the relativistic axisymmetric jets demon- 
strate the ability of our code to follow complex relativistic flows as well as the flexibility 
enough to be applied to practical astrophysical problems. These simulations show that the 
morphology and dynamics of the relativistic jets are significantly influenced by the different 
equation of state and by different compositions of a relativistic perfect gas. 

This work was supported by the National Research Foundation of Korea Grant funded 
by the Korean Government (NRF-2009-351-C00029). 
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Fig. 1. — The relativistic adiabatic index, 7 n as function of inverse temperature, £, for 
different compositions of relativistic gas. The compositions of x — (electron-positron), 0.1, 
0.3, 0.6, and 1 (electron-proton) are shown using dotted, short dashed, dot-short dashed, long 
dashed, and solid curves, respectively. The exact Synge solutions for electron-positron and 
electron-proton gases are respectively drawn as the red dotted and solid lines for comparison. 
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Fig. 2. — The quantity 7*, specific enthalpy h, and sound speed c s , as functions of inverse 
temperature £ for the different equation of state. The long dashed and short dashed lines 
correspond to the ideal gas equation of state with constant adiabatic index 7 = 5/3 and 
4/3, respectively. Results for the proposed general equation of state for electron-positron 
and electron-proton gases are shown by the dotted and solid lines, respectively. 
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Fig. 3. — (a) The two-dimensional relativistic shock tube test with the general equation of 
state for an electron-positron gas for the case of the parallel velocity in the first (moderately 
relativistic) test set. The numerical computation is performed in a square box with x = y = 
[0, 1] using 512 2 cells, and the wave structures are measured along the main diagonal line 
x = y = to 1 at time t = 0.4v^2- The numerical solutions are marked with open circles 
and the analytical solutions are plotted with solid lines, (b) Same as in (a) but for the case 
of the included tangential velocity and at time t = 0.8a/2. 
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Fig. 4. — (a) The two-dimensional relativistic shock tube test with the general equation of 
state for an electron-positron gas for the case of the parallel velocity in the second (highly 
relativistic) test set. The numerical computation is performed in a square box with x = y = 
[0, 1] using 2048 2 cells, and the wave structures are measured along the main diagonal line 
x = y = to 1 at time t = 0.4v^2- The numerical solutions are marked with open circles 
and the analytical solutions are plotted with solid lines, (b) Same as in (a) but for the case 
of the included tangential velocity and at time t = 1.8y/2. 
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Fig. 5. — (a) The two-dimensional relativistic shock reflection tests with the general equation 
of state for an electron-positron gas for the case of the parallel velocity. The numerical 
computations are performed in a square box with x = y = [0, 1] using 512 2 cells, and 
structures are measured along the main diagonal line x = y = to 1 at time t = 0.8^- The 
numerical solutions are marked with open circles and analytical solutions are plotted with 
solid lines, (b) Same as in (a) but for the case of the included tangential velocity. 
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Fig. 6. — Profiles of rest mass density p, radial velocity v r , and pressure p along the radial 
distance r in the three-dimensional relativistic blast wave tests. The numerical computations 
are performed in a cube box with x = y = z = [0, 1] using 256 3 cells, and radial structures 
are measured along the main diagonal line x — y — z — Otolat time t = 0.4. Numerical 
computations carried out with the general equation of state for an electron-proton gas and 
those done with an ideal gas equation of state with constant adiabatic index of 7 = 5/3 are 
marked using open and filled circles, respectively. 
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Fig. 7. — Images of rest mass density (top), pressure (middle), and Lorentz factor (bottom) 
in the three-dimensional relativistic blast wave test with the general equation of state for 
an electron-proton gas. The numerical computations are performed in a cube box with 
x — y — z — [0, 1] using 256 3 cells, and the images are shown in logarithmic scales on the 
plane x = at time t = 0.4. 
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Fig. 8. — Images of the rest mass density for cases A to D (top to bottom) in the three- 
dimensional simulations of relativistic axisymmetric jets. The numerical simulations have 
been performed in the computational box with x = [0,1], y = [0,1], and z = [0,4] using 
128 x 128 x 512 cells, and the images are shown in logarithmic scales on the plane x = at 
time t — 6. 
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Fig. 9. — The position of the bow shock as a function of time for cases A to D in the 
three-dimensional simulations of relativistic axisymmetric jets. The numerical results are 
marked with squares, triangles, crosses, and circles for cases A to D at selected time intervals, 
respectively, and the long dashed, short dashed, dotted, and solid lines give the corresponding 
one-dimensional theoretical estimates. 
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Table I. Norm Errors for Relativistic Shock Tube/Reflection Tests 



Test \\E(p)\\ \\E(v p )\\ \\E(vt)\\ \\E(p)\\ 



RST3a 


4.96E- 


-02 


3.91E- 


-03 


0.00E4 


-00 


3.46E- 


-02 


RST3b 


7.85E- 


-02 


2.38E- 


-03 


1.64E- 


-03 


3.11E- 


-02 


RST4a 


4.43E- 


-02 


2.83E- 


-03 


0.00E4 


-00 


5.12E- 


-01 


RST4b 


1.43E- 


-01 


2.20E- 


-03 


8.46E- 


-04 


4.21E- 


-01 


RSR5a 


1.81E- 


-01 


2.33E- 


-03 


0.00E+00 


2.61E- 


-01 


RSR5b 


2.48E- 


-01 


2.92E- 


-03 


9.92E- 


-04 


6.91E- 


-01 



Note. - The test designation indicates relativistic 
shock tube (RST) and relativistic shock reflection (RSR) 
problems followed by corresponding figure numbers (3 to 
5) and labels (a and b). 



-32 - 



Table 2. Simulation Parameters for Relativistic Axisymmetric Jets 



Case 


V 




r 6 




Pb 




EOS 


X 


iVcell 


t 


A 


10" 


-2 


7.1 


2 


2.32 x 10" 


-2 


5/3 




128 x 128 x 512 


8 


B 


10" 


-2 


7.1 


2 


6.94 x 10" 


-2 


4/3 




128 x 128 x 512 


6 


C 


10" 


-2 


7.1 


2 


3.48 x 10- 


-2 


e~e + 


0.0 


128 x 128 x 512 


7 


D 


10" 


-2 


7.1 


2 


4.23 x 10- 


-2 


e~p + 


1.0 


128 x 128 x 512 


7 



Note. — Here 77 is the density ratio of the beam to the ambient medium, 
Tb is the beam Lorentz factor, Mj, is the beam Mach number, pj, is the 
uniform pressure in both the beam and the ambient medium, EOS is the 
type of equation of state, x is the relative fraction of proton and electron 
number densities, iV ce u is the numerical resolution of the three-dimensional 
computational box, and t is the time out to which each simulation is fol- 
lowed. 



